//////////////////////////////////////////////////////////////////////
// ImgMath.h
//
// SHEN Fangyang
// me@shenfy.com
//
// Copyright (C) SHEN Fangyang, 2011, All rights reserved.
//////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////
//  Brief Description
//
//	Common math operations required in images processing.
//
//////////////////////////////////////////////////////////////////////
#pragma once
#include <cmath>

#define EPSILON1 1e-5f

namespace grflib
{
	namespace image
	{
		namespace math
		{

			inline void BilinearWeight(float x, float y, float *w)
			{
				float ww;
				float dx = x - floor(x);
				float dy = y - floor(y);
				w[0] = (1 - dx) * (1 - dy);
				w[1] = dx * (1 - dy);
				w[2] = (1 - dx) * dy;
				w[3] = dx * dy;
				ww = w[0] + w[1] + w[2] + w[3];

				for (int i = 0; i < 4; i++)
					w[i] /= ww;
			}

			inline void BicubicWeight(float x, float y, float *w)
			{
				float ww;

				float ix = x - floor(x), iy = y - floor(y);
				if(x >=0 && y >= 0 ){
					w[0] = (1-ix)*(1-iy);
					w[1] = (1-ix)*(iy);
					w[2] = (ix)*(1-iy);
					w[3] = (ix)*(iy);
					ww = w[0]+w[1]+w[2]+w[3];
				}else if( x >= 0){
					w[0] = (1-ix)*(iy+1);
					w[1] = (1-ix)*(-iy);
					w[2] = (ix)*(iy+1);
					w[3] = (ix)*(-iy);
					ww = w[0]+w[1]+w[2]+w[3];
				}else if( y >= 0){
					w[0] = (ix+1 )*(1-iy);
					w[1] = (ix+1 )*(iy);
					w[2] = (-ix)*(1-iy);
					w[3] = (-ix)*(iy);
					ww = w[0]+w[1]+w[2]+w[3];
				}else{
					w[0] = (ix+1 )*(iy+1);
					w[1] = (ix+1 )*(-iy);
					w[2] = (-ix)*(iy+1);
					w[3] = (-ix)*(-iy);
					ww = w[0]+w[1]+w[2]+w[3];
				}

				w[0] /= ww; w[1] /= ww; w[2] /= ww; w[3] /= ww;
			}

			inline void Normalize3D(float &x, float &y, float &z)
			{
				float norm = sqrt(x * x + y * y + z * z);
				if (abs(norm) > EPSILON1) //avoid numerical instability
				{
					x /= norm;
					y /= norm;
					z /= norm;
				}
			}

			__forceinline float Sign(float x)
			{
				if (x > 0.0f)
					return 1.0f;
				else if (x < 0.0f)
					return -1.0f;
				else
					return 0.0f;
			}

            ///////////////////////////////////////////////////
            //left handed coordinate system
            //y points up to the sky
			inline void SphericalToCartesianRad(float theta, float phi, float &x, float &y, float &z)
			{
                y = cos(theta);
                x = sin(theta) * cos(phi);
                z = sin(theta) * sin(phi);
			}

			inline void SphericalToCartesianDeg(float theta, float phi, float &x, float &y, float &z)
			{
				const float PI = 3.1415926535897932384f;
				theta = theta * PI / 180.0f;
				phi = phi * PI / 180.0f;
                y = cos(theta);
                x = sin(theta) * cos(phi);
                z = sin(theta) * sin(phi);
			}

			inline void CartesianToSpherical(float x, float y, float z, float &theta, float &phi)
			{
				float r = sqrt(x * x + y * y + z * z);
				theta = acos(y / r);
				phi = atan2(z, x);
			}
            ///////////////////////////////////////////////////

            //in place saturate (clamp)
            template<typename T>
            inline void SaturateIn(T &val, const T &lb, const T &ub)
            {
                if (val < lb) val = lb;
                else if (val > ub) val = ub;
            }

            template<typename T>
            inline T Saturate(const T &val, const T &lb, const T &ub)
            {
                if (val < lb)
                    return lb;
                else if (val > ub)
                    return ub;
                else
                    return val;
            }
		};
	} //namespace namespace image
} //namespace namespace grflib